Comparative analysis of sequencing technologies for single-cell transcriptomics¶
Abstract¶
All single-cell RNA-seq protocols and technologies require library preparation prior to sequencing on a platform such as Illumina. Here, we present the first report to utilize the BGISEQ-500 platform for scRNA-seq, and compare the sensitivity and accuracy to Illumina sequencing. We generate a scRNA-seq resource of 468 unique single-cells and 1,297 matched single cDNA samples, performing SMARTer and Smart-seq2 protocols on mESCs and K562 cells with RNA spike-ins. We sequence these libraries on both BGISEQ-500 and Illumina HiSeq platforms using single- and paired-end reads. The two platforms have comparable sensitivity and accuracy in terms of quantification of gene expression, and low technical variability. Our study provides a standardised scRNA-seq resource to benchmark new scRNA-seq library preparation protocols and sequencing platforms.
Data proccessing¶
SCQUA is a python package that we developed for analyzing the single cell RNA Sequencing quality.
In [1]:
%pylab inline
from SCQUA import *
In [2]:
ls ../
Figs include the resulted figures
data include all the raw data
jup include the jupyter notebook to process the data and make the figures
pdata include all the processed data.
In [3]:
ls ../data/
Salmon outputs 3 information, the counts, tpm and qc metrics. We concatenate these information in *cts.gz, *tpm.gz and *qc.gz.
detection limit (sensitivity) and accuracy based on spike-ins can be calculated.
In [3]:
ercc = get_ERCC()
sirv = get_SIRV()
spike = pd.concat([ercc,sirv])
Get results for all the datasets.
In [39]:
for f in iglob('../data/*_tpm.csv.gz'):
name = f.split('/')[-1].replace('_tpm.csv.gz',"")
print(name)
cts = pd.read_csv(f.replace('tpm','cts'), index_col=0)
tpm = pd.read_csv(f, index_col=0)
phn = pd.read_csv(f.replace('tpm','qc'), index_col=0)
df = get_result(tpm, ercc, sirv, spike)
phn = pd.concat([phn,df,cts.loc[cts.index.str.startswith("ENS")].T], axis=1)
phn["Total_counts"] = cts.loc[cts.index.str.startswith("ENS")].sum()
phn["study"] = name
phn.to_csv("../pdata/%s.csv"%name)
This file is annotation of the all the cells sequenced.
In [7]:
phn = pd.read_csv("../data/phn.csv",index_col=0)
In [3]:
phn.head()
Out[3]:
In [7]:
phn.tail()
Out[7]:
map the annotation
In [18]:
for f in iglob("../pdata/*.csv"):
print(f)
fout = f.replace("pdata","pdata1")
df = pd.read_csv(f,index_col=0)
df = df.loc[df.index.isin(phn.index)]
dfp = phn.loc[df.index]
df = pd.concat([df, dfp], axis=1)
df.to_csv(fout)
In [ ]:
Fig 1B 1C¶
(B) Single-cell detection limit (Sensitivity) of mESC cells, downsampled across two orders of magnitude from SMARTer and two Smart-seq2 replicates (633 samples). The single-cell sensitivities are largely similar between different library preparation across scRNA-seq protocols. (C) Single-cell accuracy of mESC cells, downsampled across two orders of magnitude for SMARTer and two Smart-seq2 replicates (633 samples). The grey dotted lines indicate downsampling at different read depths per cell, while red line indicates saturation per cell.
In [79]:
dfs = []
for f in iglob("../pdata1/*_ds_*"):
df = pd.read_csv(f, index_col=0)
df["id"] = df.index
df.index = df.study+"_"+df.id
df = df[["id",'Cell','Fullname','num_processed',\
'detection_limit','accuracy','Protocol','Batch','Protocol2']+
df.columns[df.columns.str.startswith('ENS')].tolist()]
dfs.append(df)
df = pd.concat(dfs,axis =0)
In [80]:
df.shape
Out[80]:
In [81]:
df.head()
Out[81]:
In [82]:
df=df[df.Protocol.str.startswith("S")]
use only the protocols in sanger
In [83]:
df = df.loc[~(df.Protocol2.str.startswith('BGISEQ-500 K562')|\
df.Protocol2.str.startswith('Hiseq4000')|\
df.Protocol2.str.startswith('BGISEQ-500 ESC'))]
In [84]:
df.Protocol.value_counts()
Out[84]:
In [85]:
df.Protocol2.value_counts()
Out[85]:
In [86]:
df=df.replace([np.inf, -np.inf], np.nan).dropna(axis=0, how='any')
In [87]:
df.to_csv("../pdata1/Fig1BC.csv")
In [2]:
df = pd.read_csv("../pdata1/Fig1BC.csv", index_col=0)
In [88]:
figsize(5,5)
matplotlib.rcParams.update({'font.size': 14})
names = ['BGISEQ-500 Smart-Seq2 rep1','HiSeq-2000 Smart-Seq2 rep1',
'BGISEQ-500 Smart-Seq2 rep2','HiSeq-2000 Smart-Seq2 rep2',
'BGISEQ-500 SMARTer','HiSeq-2000 SMARTer']
ax = fit_sensitivity(df,key2="num_processed",key3='Protocol2', \
yscale='log', xlabel = 'Sequenced reads per single-well',\
ylabel = 'Detection limit \n (# of molecules)', \
ylim = [0, 1e4], xlim = [500, 1e8], \
xlabelsize = 18, ylabelsize = 18, \
title = 'Sensitivity', titlesize = 20, names = names)
In [89]:
figsize(5,5)
matplotlib.rcParams.update({'font.size': 14})
names = ['BGISEQ-500 Smart-Seq2 rep1','HiSeq-2000 Smart-Seq2 rep1',
'BGISEQ-500 Smart-Seq2 rep2','HiSeq-2000 Smart-Seq2 rep2',
'BGISEQ-500 SMARTer','HiSeq-2000 SMARTer']
ax = fit_sensitivity(df,key2="num_processed",key3='Protocol2', \
yscale='log', xlabel = 'Sequenced reads per single-well',\
ylabel = 'Detection limit \n (# of molecules)', \
ylim = [0, 1e4], xlim = [500, 1e8], \
xlabelsize = 18, ylabelsize = 18, \
title = 'Sensitivity', titlesize = 20, names = names, colordots = True)
In [90]:
ax = fit_sensitivity(df, fun = 'accuracy',key1="accuracy",key2="num_processed",key3='Protocol2', \
xlabel = 'Sequenced reads per single-well',\
ylabel = 'Accuracy \n (pearson correlation \n coefficient)', \
ylim = [0, 1], xlim = [500, 1e8], \
xlabelsize = 18, ylabelsize = 18, \
title = 'Accuracy', titlesize = 20, names = names)
In [91]:
ax = fit_sensitivity(df, fun = 'accuracy',key1="accuracy",key2="num_processed",key3='Protocol2', \
xlabel = 'Sequenced reads per single-well',\
ylabel = 'Accuracy \n (pearson correlation \n coefficient)', \
ylim = [0, 1], xlim = [500, 1e8], \
xlabelsize = 18, ylabelsize = 18, \
title = 'Accuracy', titlesize = 20, names = names, colordots=True)
Fig 1D 1E¶
(D) PCA for matched single-cell cDNA samples performed using SMARTer and two replicates of Smart-seq2 and sequenced across both sequencing platforms. Red and green colored circles indicate sequencing of matched cDNA on Illumina HiSeq2500 and BGISEQ-500 respectively. The dotted lines represent distance i.e. measure of similarity across sequencing platforms. (E) Single-cell correlations for each scRNA-seq protocol and across sequencing platforms. The correlations (R=0.52~0.70) are comparable between sequencing platforms.
In [2]:
df = pd.read_csv("../pdata1/Fig1BC.csv", index_col=0)
In [57]:
figsize(5,5)
dfs1 = []
dfs2 = []
for p, l in zip(['Smart-Seq2 rep1','Smart-Seq2 rep2','SMARTer'],[230, 120, 210]):
df1, df2 = plot_Fig1E(df, protocol = p,\
lim = l,\
xlabel = "BGISEQ-500", ylabel="HiSeq-2000",\
xlabelsize = 18, ylabelsize=18, titlesize = 18, legendsize = 16)
plt.show()
dfs1.append(df1)
dfs2.append(df2)
[PR12&PR8]¶
In [5]:
def plot_Fig1E(df, \
protocol = "SMARTer",\
key1 = 'Batch',\
key2 = 'Protocol2',\
key3 = 'Cell',\
key4 = 'Fullname',\
lim = 210,\
cutoff = 1000, \
xlabel = None, ylabel=None,\
xlabelsize = None, ylabelsize=None, titlesize = None, legendsize = None):
dfx = df[df[key1] == protocol]
dfx = dfx[dfx.index.str.find('e6')>0]
df1 = dfx[dfx[key2].str.startswith('BGISEQ-500')].sort_values(by=key3)
df2 = dfx[dfx[key2].str.startswith('HiSeq-2000')].sort_values(by=key3)
df1 = df1[df1[key4].isin(df2[key4])]
df2 = df2[df2[key4].isin(df1[key4])]
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
dfa = pd.concat([df1,df2], axis = 0)
exprs = np.log1p(dfa.T.loc[dfa.columns.str.startswith("ENS")].astype(np.float32))
pca = PCA(n_components=50)
pca_res = pca.fit_transform(scale(exprs.T, 1))
dfa["PC1"] = pca_res[:,0]
dfa["PC2"] = pca_res[:,1]
from scipy.spatial.distance import cdist
dfb = dfa[['PC1','PC2']]
dfb1 = dfb.iloc[:int((dfb.shape[0])/2)]
dfb2 = dfb.iloc[int((dfb.shape[0])/2):]
cdst = cdist(dfb1.values,dfb2.values)
# print(np.diagonal(cdst))
xx1 = df1.detection_limit[(df2.detection_limit <cutoff).tolist()]
xx2 = df2.detection_limit[(df2.detection_limit <cutoff).tolist()]
xx3 = np.diagonal(cdst)[(df2.detection_limit <cutoff).tolist()]
slope, intercept, r_value, p_value, std_err = \
scipy.stats.linregress(xx1, xx2)
from scipy.stats import spearmanr
r_value, _ = spearmanr(xx1,xx2)
plt.plot(xx1, xx2, 'k.', label = "R=%.2f"%r_value)
plt.plot(xx1[xx3>50], xx2[xx3>50], 'r.', label='')
print("outlier cells:")
print(df2.index[np.diagonal(cdst)>50])
plt.xlim(0,lim)
plt.ylim(0,lim)
if xlabel is None: xlabel = key2
if ylabel is None: ylabel = key1
if xlabelsize is None:
plt.xlabel(xlabel)
else:
plt.xlabel(xlabel, fontsize=xlabelsize)
if ylabelsize is None:
plt.ylabel(ylabel)
else:
plt.ylabel(ylabel, fontsize=ylabelsize)
if not titlesize is None:
plt.title(protocol, fontsize = titlesize)
else:
plt.title(protocol)
if legendsize is None:
plt.legend(loc='upper right')
else:
plt.legend(loc='upper right', fontsize=legendsize)
return(df1, df2)
outlier cells are listed here:
In [6]:
figsize(5,5)
dfs1 = []
dfs2 = []
for p, l in zip(['Smart-Seq2 rep1','Smart-Seq2 rep2','SMARTer'],[230, 120, 210]):
df1, df2 = plot_Fig1E(df, protocol = p,\
lim = l,\
xlabel = "BGISEQ-500", ylabel="HiSeq-2000",\
xlabelsize = 18, ylabelsize=18, titlesize = 18, legendsize = 16)
plt.show()
dfs1.append(df1)
dfs2.append(df2)
[PR8 FigS2C]¶
In [8]:
dfs = []
for f in ['../pdata1/bgi.csv','../pdata1/sanger.csv']:
df = pd.read_csv(f, index_col=0)
df["id"] = df.index
df.index = df.study+"_"+df.id
df = df[["id",'Cell','Fullname','num_processed',\
'detection_limit','accuracy','Protocol','Batch','Protocol2']+
df.columns[df.columns.str.startswith('ENS')].tolist()]
dfs.append(df)
df = pd.concat(dfs,axis =0)
In [9]:
df=df[df.Protocol.str.startswith("S")]
In [10]:
df = df.loc[~(df.Protocol2.str.startswith('BGISEQ-500 K562')|\
df.Protocol2.str.startswith('Hiseq4000')|\
df.Protocol2.str.startswith('BGISEQ-500 ESC'))]
In [11]:
df.Protocol.value_counts()
Out[11]:
In [12]:
df.Protocol2.value_counts()
Out[12]:
In [13]:
df=df.replace([np.inf, -np.inf], np.nan).dropna(axis=0, how='any')
In [14]:
df.to_csv("../pdata1/FigS2C.csv")
In [141]:
df.shape
Out[141]:
In [ ]:
In [17]:
def plot_Fig1E(df, \
protocol = "SMARTer",\
key1 = 'Batch',\
key2 = 'Protocol2',\
key3 = 'Cell',\
key4 = 'Fullname',\
lim = 210,\
cutoff = 1000, \
xlabel = None, ylabel=None,\
xlabelsize = None, ylabelsize=None, titlesize = None, legendsize = None):
dfx = df[df[key1] == protocol]
# dfx = dfx[dfx.index.str.find('e6')>0]
df1 = dfx[dfx[key2].str.startswith('BGISEQ-500')].sort_values(by=key3)
df2 = dfx[dfx[key2].str.startswith('HiSeq-2000')].sort_values(by=key3)
df1 = df1[df1[key4].isin(df2[key4])]
df2 = df2[df2[key4].isin(df1[key4])]
xx1 = df1.detection_limit[(df2.detection_limit <cutoff).tolist()]
xx2 = df2.detection_limit[(df2.detection_limit <cutoff).tolist()]
slope, intercept, r_value, p_value, std_err = \
scipy.stats.linregress(xx1, xx2)
from scipy.stats import spearmanr
r_value, _ = spearmanr(xx1,xx2)
plt.plot(xx1, xx2, 'k.', label = "R=%.2f"%r_value)
plt.xlim(0,lim)
plt.ylim(0,lim)
if xlabel is None: xlabel = key2
if ylabel is None: ylabel = key1
if xlabelsize is None:
plt.xlabel(xlabel)
else:
plt.xlabel(xlabel, fontsize=xlabelsize)
if ylabelsize is None:
plt.ylabel(ylabel)
else:
plt.ylabel(ylabel, fontsize=ylabelsize)
if not titlesize is None:
plt.title(protocol, fontsize = titlesize)
else:
plt.title(protocol)
if legendsize is None:
plt.legend(loc='upper right')
else:
plt.legend(loc='upper right', fontsize=legendsize)
return(df1, df2)
In [18]:
figsize(5,5)
dfs1 = []
dfs2 = []
for p, l in zip(['Smart-Seq2 rep1','Smart-Seq2 rep2','SMARTer'],[230, 120, 210]):
df1, df2 = plot_Fig1E(df, protocol = p,\
lim = l,\
xlabel = "BGISEQ-500", ylabel="HiSeq-2000",\
xlabelsize = 18, ylabelsize=18, titlesize = 18, legendsize = 16)
plt.show()
dfs1.append(df1)
dfs2.append(df2)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [38]:
def plot_Fig1D(df1, df2, \
protocol = "SMARTer", \
key2 = 'Protocol2', \
key3 = 'BGISEQ-500', \
key4 = 'Fullname',\
xlabelsize = None, ylabelsize = None, titlesize = None):
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
df = pd.concat([df1,df2], axis = 0)
exprs = np.log1p(df.T.loc[df.columns.str.startswith("ENS")].astype(np.float32))
pca = PCA(n_components=50)
pca_res = pca.fit_transform(scale(exprs.T, 1))
df["PC1"] = pca_res[:,0]
df["PC2"] = pca_res[:,1]
g = df.groupby('Protocol2')
colors = ["#00cc99","#e0115f"]
for i, group in enumerate(g.indices):
tmp = g.indices[group]
plt.scatter(df.PC1.iloc[tmp], df.PC2.iloc[tmp], label=group.replace(protocol,""), lw = 0, c=colors[i])
from scipy.spatial.distance import cdist
for i in df[df[key2].str.startswith(key3)][key4]:
j = df[(df[key2].str.startswith(key3)) & (df[key4] == i)][['PC1','PC2']]
k = df[(~df[key2].str.startswith(key3)) & (df[key4] == i)][['PC1','PC2']]
# d1 = abs(j["PC1"].values - k["PC1"].values)
# d2 = abs(j["PC2"].values - k["PC2"].values)
# d = np.sqrt(d1*d1+d2*d2)
cdst = cdist(j[['PC1','PC2']].values,k[['PC1','PC2']].values)
d = cdst[0][0]
if d < 60:
plt.plot([j["PC1"],k["PC1"]],[j["PC2"],k["PC2"]],'k-.')
else:
plt.plot([j["PC1"],k["PC1"]],[j["PC2"],k["PC2"]],'r-.')
plt.xlabel("PC1 (%.3f)"%pca.explained_variance_ratio_[0], fontsize = xlabelsize)
plt.ylabel("PC2 (%.3f)"%pca.explained_variance_ratio_[1], fontsize = ylabelsize)
if not titlesize is None:
plt.title(protocol, fontsize = titlesize)
else:
plt.title(protocol)
plt.legend(scatterpoints=3, bbox_to_anchor=(1.03, 1), borderaxespad=0.);
In [39]:
i=0
protocols = ['Smart-Seq2 rep1','Smart-Seq2 rep2','SMARTer']
for df1, df2 in zip(dfs1, dfs2):
plt.clf()
plot_Fig1D(df1, df2, protocol = protocols[i], xlabelsize = 18, ylabelsize=18, titlesize = 18)
plt.show()
i+=1
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [90]:
def plot_Fig1D(df1, df2, \
protocol = "SMARTer", \
key2 = 'Protocol2', \
key3 = 'BGISEQ-500', \
key4 = 'Fullname',\
xlabelsize = None, ylabelsize = None, titlesize = None, tp = 1):
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
df = pd.concat([df1,df2], axis = 0)
exprs = np.log1p(df.T.loc[df.columns.str.startswith("ENS")].astype(np.float32))
pca = PCA(n_components=50)
pca_res = pca.fit_transform(scale(exprs.T, 1))
df["PC1"] = pca_res[:,0]
df["PC2"] = pca_res[:,1]
if tp ==1:
g = df.groupby('Protocol2')
colors = ["#00cc99","#e0115f"]
for i, group in enumerate(g.indices):
tmp = g.indices[group]
plt.scatter(df.PC1.iloc[tmp], df.PC2.iloc[tmp], label=group.replace(protocol,""), lw = 0, c=colors[i])
for i in df[df[key2].str.startswith(key3)][key4]:
j = df[(df[key2].str.startswith(key3)) & (df[key4] == i)][['PC1','PC2']]
k = df[(~df[key2].str.startswith(key3)) & (df[key4] == i)][['PC1','PC2']]
plt.plot([j["PC1"],k["PC1"]],[j["PC2"],k["PC2"]],'k-.')
plt.legend(scatterpoints=3, bbox_to_anchor=(1.03, 1), borderaxespad=0.);
else:
if tp ==2:
df1 = df[df.columns[df.columns.str.startswith('ENS')]]
df['n_genes'] = np.sum(df1>0,axis=1)
plt.scatter(df.PC1, df.PC2, lw = 0, c=df['n_genes'])
plt.colorbar()
elif tp ==3:
df['Nanog'] = df['ENSMUSG00000012396.12']
plt.scatter(df.PC1, df.PC2, lw = 0, c=df['Nanog'])
plt.colorbar()
elif tp ==4:
df['Pou5f1'] = df['ENSMUSG00000024406.16']
plt.scatter(df.PC1, df.PC2, lw = 0, c=df['Pou5f1'])
plt.colorbar()
elif tp ==5:
df['Pax6'] = df['ENSMUSG00000027168.21']
plt.scatter(df.PC1, df.PC2, lw = 0, c=df['Pax6'])
plt.colorbar()
plt.xlabel("PC1 (%.3f)"%pca.explained_variance_ratio_[0], fontsize = xlabelsize)
plt.ylabel("PC2 (%.3f)"%pca.explained_variance_ratio_[1], fontsize = ylabelsize)
if not titlesize is None:
plt.title(protocol, fontsize = titlesize)
else:
plt.title(protocol)
[PR6]¶
In [72]:
i=0
protocols = ['Smart-Seq2 rep1','Smart-Seq2 rep2','SMARTer']
for df1, df2 in zip(dfs1, dfs2):
plt.clf()
plot_Fig1D(df1, df2, protocol = protocols[i], xlabelsize = 18, ylabelsize=18, titlesize = 18,tp=2)
plt.show()
i+=1
In [73]:
df.shape
Out[73]:
In [75]:
ls ../data/GRCm38.cdna.all.symbol.tsv
In [84]:
df1 = pd.read_csv('../data/GRCm38.cdna.all.symbol.tsv', sep=' ',header=None, index_col=1)
In [86]:
df1.loc['Nanog']
Out[86]:
In [88]:
df1.loc['Pou5f1']
Out[88]:
In [89]:
df1.loc['Pax6']
Out[89]:
In [ ]:
[cells in the red circle]¶
In [70]:
df = pd.concat([dfs1[0],dfs2[0]])
In [72]:
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
exprs = np.log1p(df.T.loc[df.columns.str.startswith("ENS")].astype(np.float32))
pca = PCA(n_components=50)
pca_res = pca.fit_transform(scale(exprs.T, 1))
df["PC1"] = pca_res[:,0]
df["PC2"] = pca_res[:,1]
In [74]:
df[['PC1','PC2']].to_csv("PCA_cell_SM2rep1.csv")
In [75]:
df = pd.concat([dfs1[1],dfs2[1]])
In [76]:
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
exprs = np.log1p(df.T.loc[df.columns.str.startswith("ENS")].astype(np.float32))
pca = PCA(n_components=50)
pca_res = pca.fit_transform(scale(exprs.T, 1))
df["PC1"] = pca_res[:,0]
df["PC2"] = pca_res[:,1]
In [77]:
df[['PC1','PC2']].to_csv("PCA_cell_SM2rep2.csv")
In [ ]:
In [78]:
df = pd.concat([dfs1[2],dfs2[2]])
In [79]:
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
exprs = np.log1p(df.T.loc[df.columns.str.startswith("ENS")].astype(np.float32))
pca = PCA(n_components=50)
pca_res = pca.fit_transform(scale(exprs.T, 1))
df["PC1"] = pca_res[:,0]
df["PC2"] = pca_res[:,1]
In [80]:
df[['PC1','PC2']].to_csv("PCA_cell_SMARTer.csv")
In [ ]:
In [ ]:
In [ ]:
[PR5]¶
Nanog
In [91]:
i=0
protocols = ['Smart-Seq2 rep1','Smart-Seq2 rep2','SMARTer']
for df1, df2 in zip(dfs1, dfs2):
plt.clf()
plot_Fig1D(df1, df2, protocol = protocols[i], xlabelsize = 18, ylabelsize=18, titlesize = 18,tp=3)
plt.show()
i+=1
Pou5f1
In [92]:
i=0
protocols = ['Smart-Seq2 rep1','Smart-Seq2 rep2','SMARTer']
for df1, df2 in zip(dfs1, dfs2):
plt.clf()
plot_Fig1D(df1, df2, protocol = protocols[i], xlabelsize = 18, ylabelsize=18, titlesize = 18,tp=4)
plt.show()
i+=1
Pax6
In [93]:
i=0
protocols = ['Smart-Seq2 rep1','Smart-Seq2 rep2','SMARTer']
for df1, df2 in zip(dfs1, dfs2):
plt.clf()
plot_Fig1D(df1, df2, protocol = protocols[i], xlabelsize = 18, ylabelsize=18, titlesize = 18,tp=5)
plt.show()
i+=1
In [ ]:
In [ ]:
In [ ]:
In [ ]:
Fig 2B¶
In [31]:
dfs = []
for f in iglob("../pdata1/*_ds_*"):
df = pd.read_csv(f, index_col=0)
df["id"] = df.index
df.index = df.study+"_"+df.id
df = df[["id",'Cell','Fullname','num_processed',\
'detection_limit_ERCC','accuracy_ERCC','Protocol','Batch','Protocol2']+
df.columns[df.columns.str.startswith('ENS')].tolist()]
dfs.append(df)
df = pd.concat(dfs,axis =0)
In [32]:
df.shape
Out[32]:
In [33]:
df.head()
Out[33]:
In [34]:
df.Protocol.value_counts()
Out[34]:
In [35]:
df.Protocol2.value_counts()
Out[35]:
In [36]:
df['Protocol2'] = df['Protocol2'].str.replace("Seuqencing",'Sequencing')
In [37]:
df = df[df.Protocol2.str.find("Sequencing")>0]
In [38]:
df.Protocol2.value_counts()
Out[38]:
In [39]:
df.to_csv("../pdata1/Fig2B.csv")
In [2]:
df = pd.read_csv("../pdata1/Fig2B.csv", index_col=0)
In [3]:
df=df.replace([np.inf, -np.inf], np.nan).dropna(axis=0, how='any')
In [5]:
figsize(5,5)
matplotlib.rcParams.update({'font.size': 14})
names = ['BGISEQ-500 ESC Sequencing rep1',
'Hiseq4000 ESC Sequencing rep1',
'BGISEQ-500 K562 Sequencing rep1','Hiseq4000 K562 Sequencing rep1']
colors = [matplotlib.colors.rgb2hex(plt.get_cmap("Paired")(i)) for i in range(6,10)]
ax = fit_sensitivity(df,fun='np.log10(detection_limit_ERCC)', key1 ='detection_limit_ERCC',key2="num_processed",key3='Protocol2', \
yscale='log', xlabel = 'Sequenced reads per single-well',\
ylabel = 'Detection limit \n (# of molecules)', \
ylim = [0, 1e5], xlim = [500, 1e8], \
xlabelsize = 18, ylabelsize = 18, \
title = 'Sensitivity', titlesize = 20, names = names, colors = colors)
In [6]:
ax = fit_sensitivity(df, fun = 'accuracy_ERCC',key1="accuracy_ERCC",key2="num_processed",key3='Protocol2', \
xlabel = 'Sequenced reads per single-well',\
ylabel = 'Accuracy \n (pearson correlation \n coefficient)', \
ylim = [0, 1], xlim = [500, 1e8], \
xlabelsize = 18, ylabelsize = 18, \
title = 'Accuracy', titlesize = 20, names = names, colors = colors)
In [7]:
ax = fit_sensitivity(df, fun = 'accuracy_ERCC',key1="accuracy_ERCC",key2="num_processed",key3='Protocol2', \
xlabel = 'Sequenced reads per single-well',\
ylabel = 'Accuracy \n (pearson correlation \n coefficient)', \
ylim = [0, 1], xlim = [500, 1e8], \
xlabelsize = 18, ylabelsize = 18, \
title = 'Accuracy', titlesize = 20, names = names, colors = colors, colordots=True)
Pseudobulk¶
In [165]:
dfs = []
for f in ['bgi','sanger']:
df = pd.read_csv("../pdata1/%s.csv"%f, index_col=0)
dfs.append(df)
df = pd.concat(dfs, axis=0)
In [166]:
df.shape
Out[166]:
In [167]:
df['Protocol2'].value_counts()
Out[167]:
In [21]:
df1 = df[df.columns[df.columns.str.startswith('ENS')|df.columns.str.startswith('Protocol2')]]
df1 = df1.groupby('Protocol2').sum()
In [23]:
df1.to_csv("../pdata1/pseudobulk.csv")
In [91]:
df1 = pd.read_csv("../pdata1/pseudobulk.csv", index_col=0)
In [92]:
df1.head(2)
Out[92]:
In [170]:
df1.index
Out[170]:
In [171]:
l1 = df1.index[df1.index.str.startswith('BGI')].tolist()
l2 = df1.index[df1.index.str.startswith('Hi')]
In [54]:
l1 = ['BGISEQ-500 ESC Sequencing rep1',
'BGISEQ-500 K562 Seuqencing rep1',
'BGISEQ-500 SMARTer',
'BGISEQ-500 Smart-Seq2 rep1',
'BGISEQ-500 Smart-Seq2 rep2']
In [55]:
l2 = ['Hiseq4000 ESC Sequencing rep1',
'Hiseq4000 K562 Sequencing rep1',
'HiSeq-2000 SMARTer',
'HiSeq-2000 Smart-Seq2 rep1',
'HiSeq-2000 Smart-Seq2 rep2']
In [56]:
names = []
corrs = []
for i,j in zip(l1,l2):
cor = df1.loc[i].corr(df1.loc[j], method='pearson')
corrs.append(cor)
names.append(i.replace('BGISEQ-500 ',''))
In [57]:
df = pd.DataFrame({'name':names,'corr':corrs})
This is the pseudobulk table.
detection_limit and accuracies can only be table like this.¶
the calculation of detection_limit and accuracies need some changes in the analysis workflow.
In [58]:
df
Out[58]:
In [93]:
df = pd.DataFrame(numpy.corrcoef(df1), index = df1.index, columns = df1.index)
This is an all-to-all correlation table.
show a heatmap¶
In [94]:
df
Out[94]:
In [47]:
df.index
Out[47]:
In [185]:
cols =['BGISEQ-500 SMARTer', 'HiSeq-2000 SMARTer','BGISEQ-500 Smart-Seq2 rep1','HiSeq-2000 Smart-Seq2 rep1',
'BGISEQ-500 Smart-Seq2 rep2','HiSeq-2000 Smart-Seq2 rep2','BGISEQ-500 ESC Sequencing rep1',
'Hiseq4000 ESC Sequencing rep1', 'BGISEQ-500 K562 Seuqencing rep1',
'Hiseq4000 K562 Sequencing rep1']
In [95]:
cols =['BGISEQ-500 SMARTer', 'HiSeq-2000 SMARTer','BGISEQ-500 Smart-Seq2 rep1','HiSeq-2000 Smart-Seq2 rep1',
'BGISEQ-500 Smart-Seq2 rep2','HiSeq-2000 Smart-Seq2 rep2']
In [96]:
df = df[cols]
In [97]:
df = df.loc[cols]
[UK only plot here:]¶
In [99]:
figsize(6,6)
fig, ax = plt.subplots()
im = ax.imshow(df.values)
# We want to show all ticks...
ax.set_xticks(np.arange(df.shape[0]))
ax.set_yticks(np.arange(df.shape[0]))
# ... and label them with the respective list entries
ax.set_xticklabels(df.index)
ax.set_yticklabels(df.index)
# Rotate the tick labels and set their alignment.
plt.setp(ax.get_xticklabels(), rotation=45, ha="right",
rotation_mode="anchor")
# Loop over data dimensions and create text annotations.
for i in range(df.shape[0]):
for j in range(df.shape[0]):
text = ax.text(j, i, '%.2f'%df.iloc[i, j],
ha="center", va="center", color="w")
# ax.set_title("Harvest of local farmers (in tons/year)")
fig.tight_layout()
# plt.show()
plt.savefig("Heatmap.pdf")
In [ ]:
In [188]:
figsize(6,6)
fig, ax = plt.subplots()
im = ax.imshow(df.values)
# We want to show all ticks...
ax.set_xticks(np.arange(df.shape[0]))
ax.set_yticks(np.arange(df.shape[0]))
# ... and label them with the respective list entries
ax.set_xticklabels(df.index)
ax.set_yticklabels(df.index)
# Rotate the tick labels and set their alignment.
plt.setp(ax.get_xticklabels(), rotation=45, ha="right",
rotation_mode="anchor")
# Loop over data dimensions and create text annotations.
for i in range(df.shape[0]):
for j in range(df.shape[0]):
text = ax.text(j, i, '%.2f'%df.iloc[i, j],
ha="center", va="center", color="w")
# ax.set_title("Harvest of local farmers (in tons/year)")
fig.tight_layout()
plt.show()
In [ ]:
In [ ]:
In [ ]:
Venn plot¶
In [81]:
from matplotlib_venn import venn3
In [129]:
df = pd.read_csv("../pdata1/pseudobulk.csv", index_col=0)
In [83]:
df
Out[83]:
In [43]:
df.shape
Out[43]:
In [84]:
from matplotlib_venn import venn2
In [83]:
i=0
print(df.index[i])
venn2([set(df.loc[l1[i]][df.loc[l1[i]]>0].index), set(df.loc[l2[i]][df.loc[l2[i]]>0].index)])
Out[83]:
In [84]:
i=1
print(df.index[i])
venn2([set(df.loc[l1[i]][df.loc[l1[i]]>0].index), set(df.loc[l2[i]][df.loc[l2[i]]>0].index)])
Out[84]:
[PR4]¶
83A overlap with 84A. 83B overlap with 84B. overlap in 83 and overlap in 84.
In [85]:
# print(df.index[i])
name1 = 'BGISEQ-500 K562 Seuqencing rep1'
name2 = 'BGISEQ-500 ESC Sequencing rep1'
venn2([set(df.loc[name1][df.loc[name1]>10].index), set(df.loc[name2][df.loc[name2]>10].index)],set_labels=(name1,name2))
Out[85]:
In [87]:
# print(df.index[i])
name1 = 'Hiseq4000 K562 Sequencing rep1'
name2 = 'Hiseq4000 ESC Sequencing rep1'
venn2([set(df.loc[name1][df.loc[name1]>10].index), set(df.loc[name2][df.loc[name2]>10].index)],set_labels=(name1,name2))
Out[87]:
In [130]:
name1 = 'Hiseq4000 K562 Sequencing rep1'
name2 = 'Hiseq4000 ESC Sequencing rep1'
x1 = df.loc[name1][df.loc[name1]>10].index
x2 = df.loc[name2][df.loc[name2]>10].index
name1 = 'BGISEQ-500 K562 Seuqencing rep1'
name2 = 'BGISEQ-500 ESC Sequencing rep1'
x3 = df.loc[name1][df.loc[name1]>10].index
x4 = df.loc[name2][df.loc[name2]>10].index
In [ ]:
In [131]:
overlap = x3[x3.isin(x4)]
In [134]:
pd.Series(overlap).to_csv("overlap.csv")
In [135]:
non1 = x3[~x3.isin(x4)]
In [136]:
pd.Series(non1).to_csv("k562only_genes.csv")
In [137]:
non2 = x4[~x4.isin(x3)]
In [138]:
pd.Series(non2).to_csv("ESConly_genes.csv")
In [ ]:
In [ ]:
In [89]:
x5 = x3[x3.isin(x1)]
x6 = x4[x4.isin(x2)]
In [90]:
venn2([set(x5), set(x6)],set_labels=('K562','ESC'))
Out[90]:
In [ ]:
In [85]:
i=2
print(df.index[i])
venn2([set(df.loc[l1[i]][df.loc[l1[i]]>0].index), set(df.loc[l2[i]][df.loc[l2[i]]>0].index)])
Out[85]:
In [86]:
i=3
print(df.index[i])
venn2([set(df.loc[l1[i]][df.loc[l1[i]]>0].index), set(df.loc[l2[i]][df.loc[l2[i]]>0].index)])
Out[86]:
In [87]:
i=4
print(df.index[i])
venn2([set(df.loc[l1[i]][df.loc[l1[i]]>0].index), set(df.loc[l2[i]][df.loc[l2[i]]>0].index)])
Out[87]:
CV^2 vs. Mean¶
In [7]:
dfs = []
for f in ['bgi','sanger']:
df = pd.read_csv("../pdata1/%s.csv"%f, index_col=0)
dfs.append(df)
df = pd.concat(dfs, axis=0)
df1 = df[df.columns[df.columns.str.startswith('ENS')|df.columns.str.startswith('Protocol2')]]
In [8]:
df = df1.T
In [9]:
names = df.loc['Protocol2'].unique()
In [10]:
for s in df.loc['Protocol2'].unique():
# print(s)
df['%s_mean'%s] = np.mean(df.loc[:,df.loc['Protocol2'] == s].iloc[:-1,:],1)
df['%s_var'%s] = np.var(df.loc[:,df.loc['Protocol2'] == s].iloc[:-1,:],1)
df['%s_drop'%s] = np.sum((df.loc[:,df.loc['Protocol2'] == s].iloc[:-1,:]==0),1) \
/ df.loc[:,df.loc['Protocol2'] == s].shape[1]
df['%s_cv2'%s] = df['%s_var'%s]/(df['%s_mean'%s]**2)
df = df.iloc[:-1,:]
In [11]:
figsize(5,5)
plt.loglog()
for s in names:
plt.scatter(df['%s_mean'%s], df['%s_var'%s], s = 0.1, label = s)
plt.xlim(1e-6, 1e5)
plt.ylim(1e-9, 1e10);
plt.legend(bbox_to_anchor=(1.05, 1))
plt.xlabel("Means")
plt.ylabel("variances")
Out[11]:
[PR1] split into 5 paired ones.¶
In [12]:
names
Out[12]:
In [13]:
namet = [['BGISEQ-500 Smart-Seq2 rep2','HiSeq-2000 Smart-Seq2 rep2'],
['BGISEQ-500 Smart-Seq2 rep1','HiSeq-2000 Smart-Seq2 rep1'],
['BGISEQ-500 SMARTer','HiSeq-2000 SMARTer'],
['BGISEQ-500 K562 Seuqencing rep1','Hiseq4000 K562 Sequencing rep1',],
['BGISEQ-500 ESC Sequencing rep1','Hiseq4000 ESC Sequencing rep1']]
titles = ['Smart-Seq2 rep2',
'Smart-Seq2 rep1',
'SMARTer',
'K562',
'ESC']
In [15]:
figsize(5,5)
for tt, nameu in zip(titles,namet):
plt.clf()
plt.loglog()
for s in nameu:
plt.scatter(df['%s_mean'%s], df['%s_cv2'%s], s = 0.1, label = s)
plt.xlim(1e-6, 1e5)
plt.ylim(1e-2, 1e2);
plt.legend(bbox_to_anchor=(1.05, 1))
plt.title(tt)
plt.xlabel("Means")
plt.ylabel("$CV^2$")
plt.show()
Dropout¶
[PR2] split into 5 as well.¶
In [17]:
df.head()
Out[17]:
In [18]:
for tt, nameu in zip(titles,namet):
plt.clf()
plt.xscale('log')
for s in nameu:
plt.scatter(df['%s_mean'%s], df['%s_drop'%s], s = 0.1, label = s)
plt.xlim(1e-5, 1e5)
plt.ylim(-0.1, 1.1);
plt.legend(bbox_to_anchor=(1.05, 1))
plt.xlabel("Means")
plt.ylabel("Dropout rate")
plt.show()
Heatmap¶
In [224]:
dfs = []
for f in ['bgi','sanger']:
df = pd.read_csv("../pdata1/%s.csv"%f, index_col=0)
dfs.append(df)
df = pd.concat(dfs, axis=0)
In [225]:
df['Protocol2'] = df['Protocol2'].str.replace('Seuqencing','Sequencing')
In [102]:
dm = pd.read_csv("../data/GRCm38.cdna.all.symbol.tsv",index_col=1, sep=' ',header=None)
dn = pd.read_csv("../data/PlurinetGenes_Category.txt", sep = '\t', index_col=0)
In [103]:
dm = dm.loc[dn.index]
dm['funct_cat'] = dn['funct_cat']
dm = dm.dropna()
In [228]:
df = df[df.columns[df.columns.isin(dm[0])|df.columns.isin(['Protocol2','n_genes'])].tolist()]
In [104]:
dm['gene'] = dm.index
dm.index = dm[0]
In [105]:
dm.head()
Out[105]:
In [232]:
df.head()
Out[232]:
In [233]:
df.columns = ['n_genes']+dm.loc[df.columns[1:-1]]['gene'].tolist()+['Protocol2']
In [234]:
df.to_csv("../pdata1/Heatmap.csv")
In [106]:
dm.funct_cat.value_counts()
Out[106]:
Differentiation¶
In [100]:
df = pd.read_csv("../pdata1/Heatmap.csv", index_col=0)
In [107]:
dm[dm.funct_cat == 'Differentiation']['gene'].tolist() + ['Protocol2','n_genes']
Out[107]:
In [237]:
df.Protocol2.value_counts()
Out[237]:
In [108]:
df = df[dm[dm.funct_cat == 'Differentiation']['gene'].tolist() + ['n_genes','Protocol2']]
left side is BGI
In [127]:
def make_heatmap(df, key = 'HiSeq-2000 SMARTer', cut=5000,title=''):
x = df[df.Protocol2.str.endswith(key)]
z = x.iloc[:,:-1]
z = z.sort_values(by='n_genes')
y = z.iloc[:,:-1]
plt.clf()
plt.pcolormesh(np.log1p(y.T),cmap='RdYlBu_r')
plt.ylabel('Genes', fontsize= 24)
plt.xlabel('Cells', fontsize= 24)
plt.axvline((x['n_genes']<cut).sum(), lw=2, c='r',linestyle='-')
plt.gca().set_yticks(np.arange(y.shape[1])+0.5, minor=True)
plt.gca().set_yticks([], minor=False)
plt.gca().set_yticklabels(x.columns[:-2], minor=True)
plt.gca().set_xticks([], minor=False)
plt.title(title)
plt.show()
red line is the cutoff of number of genes. 7500 for BGI and 5000 for Uk.
In [122]:
x = df[df.Protocol2.str.endswith('BGISEQ-500 ESC Sequencing rep1')]
In [123]:
z = x.sort_values(by='n_genes')
In [ ]:
In [128]:
figsize(5,3)
for k in df.Protocol2.unique():
cut = 5000
if k.find('BGI')>-1:cut=7500
make_heatmap(df, key = k, cut=cut,title=k)
In [ ]:
Pluripotency¶
In [270]:
df = pd.read_csv("../pdata1/Heatmap.csv", index_col=0)
In [271]:
dm.funct_cat.value_counts()
Out[271]:
In [272]:
df.Protocol2.value_counts()
Out[272]:
In [273]:
df = df[dm[dm.funct_cat == 'Pluripotency']['gene'].tolist() + ['n_genes','Protocol2']]
In [275]:
figsize(5,5)
for k in df.Protocol2.unique():
cut = 5000
if k.find('BGI')>0:cut=7500
make_heatmap(df, key = k, cut=cut,title=k)
In [ ]:
In [240]:
figsize(5,5)
make_heatmap(df, key = 'Smart-Seq2 rep1')
In [241]:
figsize(5,5)
make_heatmap(df, key = 'Smart-Seq2 rep2')
In [242]:
figsize(5,5)
make_heatmap(df, key = 'SMARTer')
In [243]:
figsize(5,5)
make_heatmap(df, key = 'K562 Sequencing rep1')
In [244]:
figsize(5,5)
make_heatmap(df, key = 'ESC Sequencing rep1')
Non-overlapping genes¶
In [34]:
dfs = []
for f in ['bgi','sanger']:
df = pd.read_csv("../pdata1/%s.csv"%f, index_col=0)
dfs.append(df)
df = pd.concat(dfs, axis=0)
In [35]:
df['Protocol2'] = df['Protocol2'].str.replace('Seuqencing','Sequencing')
In [46]:
df1 = pd.read_csv("../pdata1/pseudobulk.csv", index_col=0)
In [36]:
key = 'Smart-Seq2 rep1'
df2 = df1[df1.index.str.endswith(key)]
gl = df2.columns[((df2.iloc[0]>0)&(df2.iloc[1]==0))].tolist()+ \
df2.columns[((df2.iloc[1]>0)&(df2.iloc[0]==0))].tolist()
df3 = df[gl + ['Protocol2']]
df2 = df3[df3.Protocol2.str.endswith(key)]
df2 = df2.iloc[:,list(arange(100))+list(arange(-1, -102, -1))[::-1]]
figsize(10,10)
make_heatmap(df2, key = key)
In [37]:
key = 'Smart-Seq2 rep2'
df2 = df1[df1.index.str.endswith(key)]
gl = df2.columns[((df2.iloc[0]>0)&(df2.iloc[1]==0))].tolist()+ \
df2.columns[((df2.iloc[1]>0)&(df2.iloc[0]==0))].tolist()
df3 = df[gl + ['Protocol2']]
df2 = df3[df3.Protocol2.str.endswith(key)]
df2 = df2.iloc[:,list(arange(100))+list(arange(-1, -102, -1))[::-1]]
figsize(10,10)
make_heatmap(df2, key = key)
In [38]:
key = 'SMARTer'
df2 = df1[df1.index.str.endswith(key)]
gl = df2.columns[((df2.iloc[0]>0)&(df2.iloc[1]==0))].tolist()+ \
df2.columns[((df2.iloc[1]>0)&(df2.iloc[0]==0))].tolist()
df3 = df[gl + ['Protocol2']]
df2 = df3[df3.Protocol2.str.endswith(key)]
df2 = df2.iloc[:,list(arange(100))+list(arange(-1, -102, -1))[::-1]]
figsize(10,10)
make_heatmap(df2, key = key)
In [47]:
df1.index = df1.index.str.replace('Seuqencing','Sequencing')
In [48]:
key = 'K562 Sequencing rep1'
df2 = df1[df1.index.str.endswith(key)]
gl = df2.columns[((df2.iloc[0]>0)&(df2.iloc[1]==0))].tolist()+ \
df2.columns[((df2.iloc[1]>0)&(df2.iloc[0]==0))].tolist()
df3 = df[gl + ['Protocol2']]
df2 = df3[df3.Protocol2.str.endswith(key)]
df2 = df2.iloc[:,list(arange(100))+list(arange(-1, -102, -1))[::-1]]
figsize(10,10)
make_heatmap(df2, key = key)
In [41]:
key = 'ESC Sequencing rep1'
df2 = df1[df1.index.str.endswith(key)]
gl = df2.columns[((df2.iloc[0]>0)&(df2.iloc[1]==0))].tolist()+ \
df2.columns[((df2.iloc[1]>0)&(df2.iloc[0]==0))].tolist()
df3 = df[gl + ['Protocol2']]
df2 = df3[df3.Protocol2.str.endswith(key)]
df2 = df2.iloc[:,list(arange(100))+list(arange(-1, -102, -1))[::-1]]
figsize(10,10)
make_heatmap(df2, key = key)
Insert size¶
In [2]:
df1 = pd.read_csv("../data/fLen_bgi.csv", index_col=0)
df2 = pd.read_csv("../data/fLen_sanger.csv", index_col=0)
df = pd.concat([df1,df2])
df.index = df.index.str.split('/').str[-1]
In [3]:
df.head()
Out[3]:
In [8]:
df['Protocol2'] = phn.loc[df.index]['Protocol2']
df = df[~df.Protocol2.isna()]
In [9]:
df.shape
Out[9]:
In [10]:
df.columns = list(arange(1001)+1)+['Protocol2']
In [11]:
df['Protocol2'] = df['Protocol2'].str.replace('Seuqencing','Sequencing')
In [12]:
df.Protocol2.value_counts()
Out[12]:
In [13]:
colors = [matplotlib.colors.rgb2hex(plt.get_cmap("Paired")(i)) for i in range(0,10)]
In [14]:
cc = dict(zip(df.Protocol2.unique(),colors))
In [17]:
figsize(5,5)
for i in arange(df.shape[0]):
plt.scatter(arange(1001)+1, df.iloc[i][:-1], s = 0.1, label = df.iloc[i][-1], c = cc[df.iloc[i][-1]])
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
plt.xticks([0,200,400,600,800,1000],[0,20,40,60,80,100])
Out[17]:
In [ ]:
In [18]:
figsize(5,5)
df1 = df[df.Protocol2.str.endswith('Smart-Seq2 rep1')]
for i in arange(df1.shape[0]):
plt.scatter(arange(1001)+1, df1.iloc[i][:-1], s = 0.1, label = df1.iloc[i][-1], c = cc[df1.iloc[i][-1]])
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
plt.xticks([0,200,400,600,800,1000],[0,20,40,60,80,100])
Out[18]:
In [19]:
figsize(5,5)
df1 = df[df.Protocol2.str.endswith('Smart-Seq2 rep2')]
for i in arange(df1.shape[0]):
plt.scatter(arange(1001)+1, df1.iloc[i][:-1], s = 0.1, label = df1.iloc[i][-1], c = cc[df1.iloc[i][-1]])
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
plt.xticks([0,200,400,600,800,1000],[0,20,40,60,80,100])
Out[19]:
In [20]:
figsize(5,5)
df1 = df[df.Protocol2.str.endswith('SMARTer')]
for i in arange(df1.shape[0]):
plt.scatter(arange(1001)+1, df1.iloc[i][:-1], s = 0.1, label = df1.iloc[i][-1], c = cc[df1.iloc[i][-1]])
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
plt.xticks([0,200,400,600,800,1000],[0,20,40,60,80,100])
Out[20]:
In [21]:
figsize(5,5)
df1 = df[df.Protocol2.str.endswith('ESC Sequencing rep1')]
for i in arange(df1.shape[0]):
plt.scatter(arange(1001)+1, df1.iloc[i][:-1], s = 0.1, label = df1.iloc[i][-1], c = cc[df1.iloc[i][-1]])
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
plt.xticks([0,200,400,600,800,1000],[0,20,40,60,80,100])
Out[21]:
In [22]:
figsize(5,5)
df1 = df[df.Protocol2.str.endswith('K562 Sequencing rep1')]
for i in arange(df1.shape[0]):
plt.scatter(arange(1001)+1, df1.iloc[i][:-1], s = 0.1, label = df1.iloc[i][-1], c = cc[df1.iloc[i][-1]])
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
plt.xticks([0,200,400,600,800,1000],[0,20,40,60,80,100])
Out[22]:
Coverage¶
In [23]:
df1 = pd.read_csv("../data/coverage_bgi.csv", index_col=0)
df2 = pd.read_csv("../data/coverage_sanger.csv", index_col=0)
phn = pd.read_csv("../data/phn.csv",index_col=0)
df = pd.concat([df1,df2], axis=0)
df.index = df.index.str.split('/').str[-1]
In [24]:
df.head()
Out[24]:
In [25]:
df['Protocol2'] = phn.loc[df.index]['Protocol2']
df = df[~df.Protocol2.isna()]
In [26]:
df.shape
Out[26]:
In [27]:
df.columns = list(np.linspace(0, 100, 20))+['Protocol2']
In [28]:
df['Protocol2'] = df['Protocol2'].str.replace('Seuqencing','Sequencing')
In [29]:
df.Protocol2.value_counts()
Out[29]:
In [30]:
colors = [matplotlib.colors.rgb2hex(plt.get_cmap("Paired")(i)) for i in range(0,10)]
In [31]:
cc = dict(zip(df.Protocol2.unique(),colors))
In [33]:
figsize(5,5)
for i in arange(df.shape[0]):
plt.plot(np.linspace(0, 100, 20), df.iloc[i][:-1], label = df.iloc[i][-1], c = cc[df.iloc[i][-1]])
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
# plt.xticks([0,200,400,600,800,1000],[0,20,40,60,80,100])
Out[33]:
[PR10]¶
In [34]:
figsize(5,5)
df1 = df[df.Protocol2.str.endswith('Smart-Seq2 rep1')]
for i in arange(df1.shape[0]):
plt.plot(np.linspace(0, 100, 20), df1.iloc[i][:-1], label = df1.iloc[i][-1], c = cc[df1.iloc[i][-1]],\
alpha=0.2)
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
Out[34]:
In [35]:
figsize(5,5)
df1 = df[df.Protocol2.str.endswith('Smart-Seq2 rep2')]
for i in arange(df1.shape[0]):
plt.plot(np.linspace(0, 100, 20), df1.iloc[i][:-1], label = df1.iloc[i][-1], c = cc[df1.iloc[i][-1]],\
alpha=0.2)
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
Out[35]:
In [36]:
figsize(5,5)
df1 = df[df.Protocol2.str.endswith('SMARTer')]
for i in arange(df1.shape[0]):
plt.plot(np.linspace(0, 100, 20), df1.iloc[i][:-1], label = df1.iloc[i][-1], c = cc[df1.iloc[i][-1]],\
alpha=0.2)
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
Out[36]:
In [37]:
figsize(5,5)
df1 = df[df.Protocol2.str.endswith('ESC Sequencing rep1')]
for i in arange(df1.shape[0]):
plt.plot(np.linspace(0, 100, 20), df1.iloc[i][:-1], label = df1.iloc[i][-1], c = cc[df1.iloc[i][-1]],\
alpha=0.2)
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
Out[37]:
In [38]:
figsize(5,5)
df1 = df[df.Protocol2.str.endswith('K562 Sequencing rep1')]
for i in arange(df1.shape[0]):
plt.plot(np.linspace(0, 100, 20), df1.iloc[i][:-1], label = df1.iloc[i][-1], c = cc[df1.iloc[i][-1]],\
alpha=0.2)
plt.xlabel('length %', fontsize=18)
plt.ylabel('density', fontsize=18)
Out[38]:
In [ ]:
Alternative Splicing Event¶
In [139]:
phn = pd.read_csv("../data/phn.csv",index_col=0)
df = pd.read_csv("../data/ASE_stat.csv", index_col=0, sep='\t')
df.columns = ['ASE']
df = df.loc[df.index.isin(phn.index)]
phn['ASE'] = 0
ds = phn['ASE']
ds.loc[df.index] = df['ASE'].tolist()
phn['ASE'] = ds
In [140]:
phn['Protocol2'] = phn['Protocol2'].str.replace("Seuqencing",'Sequencing')
In [141]:
phn.Protocol2.value_counts()
Out[141]:
In [142]:
df = phn[(phn['Protocol2'].str.startswith('BGISEQ-500 K562')|phn['Protocol2'].str.startswith('Hiseq4000 K562'))]
In [143]:
df = df[df['Protocol2'].str.endswith('rep1')]
In [144]:
df.Protocol2.value_counts()
Out[144]:
In [ ]:
In [145]:
ax=sns.violinplot(x='Protocol2',y='ASE',data=df,inner=None, color=".9",cut=0)
ax.set_xticklabels(ax.get_xticklabels(),rotation=90)
ax = sns.stripplot(x='Protocol2',y='ASE',data=df, jitter=0.3)
In [147]:
df['Read_type'].value_counts()
Out[147]:
In [148]:
ax=sns.violinplot(x='Read_type',y='ASE',data=df,inner=None, color=".9",cut=0)
ax.set_xticklabels(ax.get_xticklabels(),rotation=90)
ax = sns.stripplot(x='Read_type',y='ASE',data=df, jitter=0.3)
e6
In [13]:
phn = pd.read_csv("../data/phn.csv",index_col=0)
df = pd.read_table("../data/AS_stat.xls", index_col=0, sep='\t')
df.columns = ['ASE']
In [14]:
df = df.loc[df.index.isin(phn.index)]
phn['ASE'] = 0
ds = phn['ASE']
ds.loc[df.index] = df['ASE'].tolist()
phn['ASE'] = ds
In [15]:
ds.sum()
Out[15]:
In [ ]:
In [3]:
ls ../data/AS*
In [ ]:
In [16]:
phn['Protocol2'] = phn['Protocol2'].str.replace("Seuqencing",'Sequencing')
In [17]:
phn.Protocol2.value_counts()
Out[17]:
In [9]:
# df = phn[(phn['Protocol2'].str.startswith('BGISEQ-500 K562')|phn['Protocol2'].str.startswith('Hiseq4000 K562'))]
In [20]:
df = phn
In [21]:
df = df[df['Protocol2'].str.endswith('rep1')]
In [11]:
df.Protocol2.value_counts()
Out[11]:
In [22]:
ax=sns.violinplot(x='Read_type',y='ASE',data=df,inner=None, color=".9",cut=0)
ax.set_xticklabels(ax.get_xticklabels(),rotation=90)
ax = sns.stripplot(x='Read_type',y='ASE',data=df, jitter=0.3)
In [ ]:
In [36]:
df = phn[(phn['Protocol2'].str.startswith('BGISEQ-500 ESC')|phn['Protocol2'].str.startswith('Hiseq4000 ESC'))]
In [37]:
df = df[df['Protocol2'].str.endswith('rep1')]
In [38]:
df.Protocol2.value_counts()
Out[38]:
In [39]:
ax=sns.violinplot(x='Protocol2',y='ASE',data=df,inner=None, color=".9",cut=0)
ax.set_xticklabels(ax.get_xticklabels(),rotation=90)
ax = sns.stripplot(x='Protocol2',y='ASE',data=df, jitter=0.3)
In [22]:
phn = pd.read_csv("../pdata/bgi_human_phn.csv", index_col=0)
In [33]:
phn = pd.read_csv("../pdata/bgi_mouse_phn.csv", index_col=0)
In [34]:
phn.head()
Out[34]:
In [35]:
phn.Protocol2.value_counts()
Out[35]:
In [36]:
phn = phn[phn.Protocol2.str.startswith('BGISEQ-500 K562')|phn.Protocol2.str.startswith('BGISEQ-500 ESC')|\
phn.Protocol2.str.startswith('Hiseq4000 K562')|phn.Protocol2.str.startswith('Hiseq4000 ESC')]
In [37]:
phn.shape
Out[37]:
In [38]:
phn.Protocol2.value_counts()
Out[38]:
In [20]:
ls ../pdata
In [21]:
ls ../pdata1
In [39]:
df = pd.read_csv("../pdata/bgi.csv", index_col=0)
In [40]:
df.head()
Out[40]:
In [41]:
df.shape
Out[41]:
In [ ]:
In [42]:
df = df.loc[phn.index]
In [43]:
df.shape
Out[43]:
In [44]:
df = df[~df.num_processed.isna()]
In [45]:
df
Out[45]:
In [46]:
df1 = df.replace([np.inf, -np.inf], np.nan).dropna(subset=["detection_limit_SIRV"], how="all")
In [47]:
df1
Out[47]:
In [ ]:
In [2]:
ls ../supinfo/
In [3]:
from glob import iglob
In [7]:
for f in iglob('../supinfo/Tab*.csv'):
df = pd.read_csv(f,index_col=0)
fo = f.replace(".csv",'.1.csv')
df = df[df.columns[~df.columns.str.startswith('ENS')]]
df.to_csv(fo)
In [ ]:
In [ ]:
End¶
In [ ]: